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ABSTRACT 

We describe similarity solutions that characterize the collapse of collisional gas 
onto scalefree perturbations in an Einstein-de Sitter universe. We consider the effects 
of radiative cooling and derive self-similar solutions under the assumption that the 
cooling function is a power law of density and temperature, A(T, p) oc p^^^T. We use 
these results to test the ability of Smooth Particle Hydrodynamics (SPH) techniques 
to follow the collapse and accretion of shocked, rapidly cooling gas in a cosmological 
context. Our SPH code reproduces the analytical results very well in cases that include 
or exclude radiative cooling. No substantial deviations from the predicted central mass 
accretion rates or from the temperature, density, and velocity profiles are observed in 
well resolved regions inside the shock radius. This test problem lends support to the 
reliability of SPH techniques to model the complex process of galaxy formation. 



1 INTRODUCTION 



Structure forms in hierarchically clustering universes as pri- 
mordial dark matter density fluctuations are amplified by 
gravity and collapse in a constantly evolving population of 
virialized dark matter halos. In this scenario galaxies are en- 
visioned to form as baryons follow the dark matter collapse, 
dissipate their kinetic energy through shocks, and radiate 
it away as they settle (and form stars) in centrifugally sup- 
ported structures at the center of dark halos. Galaxies evolve 
afterwards as a result of mergers between protogalaxies 
and of further accretion of intergalactic gas (White & Rees 
1978, Navarro & White 1993, Gen & Ostriker 1993, Navarro, 
Frenk & White 1994, Evrard, Summers & Davis 1994, Katz, 
Weinberg & Hernquist 1996, Bryan et al 1998, Gouchman, 
Thomas, & Pearce 1995, Yepes et al 1997, Navarro & Stein- 
metz 1997, Tissera, Lambas & Abadi 1997, Steinmetz & 
Navarro 1998). 

Gravity, pressure gradients, hydrodynamical shocks, 
and the ability of gas to radiate are therefore physical pro- 
cesses of crucial importance during the formation of galaxies 
in a cosmological context. Numerical experiments intended 
to simulate galaxy formation must therefore capture accu- 
rately these essential ingredients on the many different levels 
of the hierarchy that coexist at a given time. Unfortunately, 
detailed analytic solutions are not known for relevant ana- 
logues of the complex galaxy formation process and it has 
been difficult to assess properly the accuracy and reliability 
of these codes (Frenk et al 1999). 

Previous studies have therefore focussed on the sen- 
sitivity and convergence of the results regarding numeri- 
cal parameters such as the size of the grid used in Eule- 



rian hydrodynamical methods (Gen 1992, Bryan et al 1998) 
or the number of particles used in particle-based methods 
such as the Smooth Particle Hydrodynamics (SPH, see Gin- 
gold & Monaghan 1977, Benz 1990, Hernquist & Katz 1989, 
Navarro & White 1993, Summers 1993 for general intro- 
ductions to SPH). Convergence as resolution improves is a 
necessary condition for simulations to be reliable, but is of- 
ten not sufficient to ensure that the results are accurate and 
robust. There is clearly a need for analytic solutions that 
describe physical situations similar to the galaxy formation 
scenario envisioned in cosmological models and that can be 
used to gauge the performance of cosmological hydrodynam- 
ical codes. 

Spherical infall is one relevant situation for which a de- 
tailed solution is known. Bertschinger (1985) first computed 
the detailed behaviour of collisional gas being accreted onto 
a point mass perturber in an Einstein-de Sitter universe. 
Assuming that only gravity, pressure gradients, and hydro- 
dynamical shocks control the gas behaviour, Bertschinger 
exploited the scalefree nature of all these processes to derive 
similarity solutions that offer a useful testbed for hydrody- 
namical codes (Navarro & White 1993, Summers 1993). One 
crucial ingredient of the galaxy formation process is, how- 
ever, missing from these tests: radiative cooling. This is be- 
cause the cooling function is the result of atomic processes 
that are not independent of scale and therefore similarity 
solutions of the spherical infall problem are not generally 
available when radiative energy losses are included. This is 
true even in the very simplified case when radiation transfer, 
heat conduction, and magnetic effects are neglected. 

Self-similarity may still be recovered in situations that 
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include radiative cooling at the expense of placing restric- 
tions on the temperature and density dependence of the 
cooling function. For example, Bertschinger (1989) com- 
puted the detailed self-similar evolution of cooling flows in 
isothermal potentials under the assumption that the gas 
cooling function is a power law of density and temperature, 
A(p,T) oc p^T^. For A = 1/2, this power law resembles 
the contribution from thermal bremsstrahlung to the over- 
all cooling function, so these results can be usefully applied 
to the hot, diffuse X-ray emitting gas that fills the intra- 
cluster medium of rich galgixy clusters. Unfortunately, these 
solutions are only valid for isolated systems originally in 
hydrostatic equilibrium and therefore their applicability to 
problems where contirmous mass accretion play a significant 
role is limited. 

A related approach has been recently described by 
Owen, Weinberg &: Villumsen (1988), who derive a fam- 
ily of cooling functions that ensure self-similar evolution in 
Einstein-de Sitter universes with power-law initial density 
fluctuations. In this case, similarity is preserved by ensuring 
that the cooling timescale of an object with characteristic 
clustering mass (M*) is a flxed fraction of the Hubble time. 

In this paper we follow a similar approach and derive 
similarity solutions for the spherical infall problem that in- 
clude energy losses due to radiative cooling. Similarity is 
preserved by choosing a convenient power-law form of the 
cooling function, and its solutions are compared with the re- 
sults of direct numerical simulations using a hydrodynamical 
SPH code. Because of the restrictions placed on the cooling 
function, the applicability of the results to realistic mod- 
els of galaxy formation is not straightforward, but the solu- 
tions are very useful as tests of hydrodynamical codes un- 
der physical conditions that combine the major ingredients 
of galaxy formation models: gravitational collapse, pressure 
gradients, energy dissipation through shocks, and radiative 
energy losses, in a proper cosmological context. 

We derive the similarity solutions in §2 and compare 
them with numerical simulations in §3. Section 4 discusses 
the results and §5 summarizes our main conclusions. 



increases with time as 

rtaoct^ (2) 

where ^ = (2/3) (1 + l/3e). The mass inside the turnaround 
radius then grows as Mta = M[r < rta) oc t'^''^'^ cx (1 -|- 
z)~^^'^, where z is the usual definition of redshift. 

This scaling can be compared with that of the charac- 
teristic clustering mass in a scalefree hierarchical clustering 
universe where the power spectrum of initial density fluc- 
tuations is P{k) oc k": AL{z) oc (1 + z)-^^^"+^\ Perturba- 
tions characterized by a given value of e therefore accrete 
mass at the same rate as a "typical" mass concentration 
in a scalefree universe with n = 3(2e — 1). A central point 
mass perturbation corresponds to e = 1, with rta oc t^^^ 
and Mta oc t^^^ oc {1 + z)~^. This is the case considered by 
Bertschinger (1985). 

Assuming that at t = ti the initial velocity fleld is pure 
unperturbed Hubble flow, v — Hir = {2/3)r/ti, there are 
no further scales in the problem once the magnitude of ini- 
tial density perturbation has been specified, and the time 
evolution of the system must approach self-similarity after 
a short initial transient. This implies that a unique solu- 
tion, oxprossod in properly scaled variables, describes the 
structure of system at all times t S> ii- It is convenient 
to express this solution in nondimensional form and, fol- 
lowing Bertschinger (1985), we define dimensionless radii, 
velocities, densities, pressures, masses, and temperatures as 
follows, 

X{r,t) = — (3.1) 

rta 

v{r,t)='-fV{X) (3.2) 
p{r,t) = pH{t)D{X) (3.3) 

p{r,t) = pHit)(^^y P{\) (3.4) 
m{r,t)^^PHrlM{X) (3.5) 



2 ACCRETION OF COLLISIONAL GAS ONTO 
SCALE FREE PERTURBATIONS 

2.1 Similarity solutions neglecting radiative 
cooling 

The evolution of a density perturbation in an Einstein-de 
Sitter universe is expected to be self-similar if the initial 
perturbation contains no physical scales. Following Fillmore 
& Goldrcich (1984) wo characterize the initial perturbation 
at some initial time ti by the excess mass within a radial 
shell of (initial) radius r;. 



6Mi 



\MoJ 



(1) 



where M, = (4/3)7rpjfrf is the unperturbed mass within 
ri, pH = 3Hf /8itG is the critical density for closure. Hi 
is Hubble's constant a,t t — ti, Mo is some reference mass, 
and e > 0. Because the mass excess is positive each radial 
shell is bound to the center and collapses after reaching a 
(maximum) turnaround radius, rta- The turnaround radius 



T(A) 



P{X) 



(3.6) 



(7-l)0(A)' 

where 7 = 5/3 is the usual ratio of specific heats. 

The equations describing the motion of a coUisional 
fluid with spherical symmetry can be expressed in terms 
of these dimensionless variables and are given by (see 
Bertschinger 1985 for details). 



{V- 



^X)D' + DV' -h -2D = Q 



{V-iX)V'-{l 



^' D 9A2 



(V 



2(2 - - 27 



M' 



SX^D 



(4.1) 



(4.2) 



(4.3) 



(4.4) 



Here primes refer to differentiation relative to A. Eqs.(4) 
are, respectively, the continuity, Euler, adiabatic, and mass 
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equations and are valid for a pressurized fluid flow neglecting 
radiation, heat conduction, and deviations from spherical 
symmetry. 

Pressurization of each radial shell of fluid occurs soon 
after turnaround as the collapsing shell encounters previ- 
ously collapsed ones. Because of similarity constraints, the 
radius at which the shock occurs must be a constant fraction 
of the turnaround radius, A = As- Outside As the evolution 
of the gas is identical to the turnaround and collapse of a 
pressureless shell of material. A full solution of eqs.(4) can 
bo found by locating the radius of the shock and applying 
Hugoniot shock jump conditions to the exterior pressure- 
less infall values. A parametric form of the preshock cold 
accretion flow is given by. 



\ = 810^(6/2) ( j 



V{\) = A 



sin 0{0 - ahiO 
(1- cos 61)2 



(5.1) 



(5.2) 



D{X) = 



{9- 



2(l-cos6')3(l + 3ex) 



3 9^^^sin^ 
^^^)-^ 2(l-cos^)3 



(5.3) 



(5.4) 



withx = l-(3/2)(l^(A)/A). 

The shock location depends also on the central bound- 
ary conditions, which wo take to be that the velocity and 
mass must vanish, i.e. y = M = OatA = 0. Identifying 
the values of the variables inside (outside) the shock by the 
subscript 2 (1), we have, at A = As, 



7- 



{V2 - eAs) 

(V^i-^As) 7+1 



D2 



7-1 
7 + 1 



P2 = 



7+1 



(6.1) 
(6.2) 
(6.3) 
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Figure 2. Examples of the three types of similarity solutions found when radiative cooling is included, computed for the case e = 1 
and Ka = 0.1. The stagnation, adiabatic, and eigensolution correspond to the dotted (As = 0.23), dashed (As = 0.15), and solid line 
(As = 0.1869), respectively. The eigensolution represents a limiting case of the two other kinds, when the stagnation point approaches 
the center and the flow extends all the way to A = 0. 



M2 = Ml. 



(6.4) 



Figure 1 shows the resulting density, velocity, temper- 
ature and entropy profiles for various values of tiie initial 
perturbation parameter e. As e decreases from unity (the 
value corresponding to a point mass perturbation, see solid 
line) to 1/3 (dashed line), the shock moves inwards, the in- 
ner density profile becomes shallower and the temperature 
profile becomes approximately isothermal. We shall see next 
how these results are altered by the inclusion of radiative 
cooling eflfects. 



2.2 Similarity solutions including radiative energy 
losses 

2.2.1 The self-similar cooling function 

The results discussed in the previous subsection are only 
applicable in the limiting case when energy losses due to 
radiative cooling are neglected. As discussed in §1, the cool- 



ing function of a plasma with realistic cosmic abundances 
has a complex dependence on temperature and imposes di- 
mensional physical scales on the problem that violate the 
conditions required for the existence of self-similar solutions. 

Similarity solutions may exist only when the cooling 
processes introduce no further scales in the problem. This 
condition can be satisfied by choosing an appropriate cooling 
function so that the overall cooling efficiency is independent 
of time. This can be ensured by demanding, for example, 
that the cooling radius (i.e. the radius at which the cooling 
time equals the age of the universe) be a fixed fraction of 
the turnaround radius of the system. Equivalently, one may 
require that, at some fixed fraction of the turnaround radius, 
the ratio between the local cooling timescale and the age of 
the universe be constant and independent of time. 

The cooling time is given by 



icool — 



up 



UpH 



(7) 



du/dt A(p,T) A(p,T)' 
where w oc T is the specific thermal energy of the gas, and 
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Figure 3. The density, veloeity, temperature, and enclosed mass profiles showing the effeet of increasing the relative importance of 
cooling. All curves correspond to the "eigensolution" for e = 1. Shock radii are given in Table 1. The solid, dotted, and dashed curves 
correspond to Kq = 0.0, 0.1 and 0.3, respectively. 



the proportionality in eq.(7) is valid at a fixed value of A. 

The condition 

tcooi/tH = {6TTGpHy^^tcooi = Kq^ = Constant (8) 

is thus satisfied if 

A(p, T) oc p'/'m cx p^/^T. (9) 

This condition is independent of e and implies that the so- 
lution will be self similar regardless of the time dependence 
of the turnaround radius. 

Our similarity solutions thus require a weaker depen- 
dence on density and a stronger dependence on tempera- 
ture than expected from thermal bremsstrahlung emission, 
^bremss OC p^T^^^ . We uote, however, that eq.(9) is not the 
only cooling function that would lead to self-similar evolu- 
tion. In particular, since the characteristic "virial temper- 
ature" {Tyir cx GMta/rta) of thc systcm is related to its 
mean density by the growth rate of the turnaround radius, 
it is possible to retain the dependence characteristic of 
realistic cooling functions and adjust only the temperature 



exponent to preserve similarity. The price one pays is that in 
this case the temperature exponent of the self-similar cool- 
ing function depends on e. More explicitly, A(p, T) oc p^T^ , 
with/3 = l-(9/2)(e/(3e-2)) (Owen et al 1998). In this case, 
the relative velocities of the cooling radius and the shock ra- 
dius are equal. We emphasize that our choice of self-similar 
cooling function (eq.9) is independent of t and does not rely 
on tuning the two velocities to agree with each other. 



2.2.2 The similarity solutions 

Once the appropriate form of the cooling function has been 
chosen, the behaviour of the gas can be computed using 
eqs.(4) after modifying the entropy conservation equation 
(4.3) to allow for energy losses. In dimensionless form, the 
modified eq.(4.3) now reads 

{V - CA) - 7^] = 2(2 - - 27 - KoD''\ (10) 
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Figure 4. Comparison between similarity solutions and the results of SPH simulations. The solution without cooling is represented by 
the solid line. We also show, for comparison, the solution including radiative energy losses [Kq = 0.1). Different symbols correspond to 
the SPH simulation at different times. Open triangles, squares and filled circles correspond to times when ~ 7, 15, and 20% of the initial 
mass lies inside the shock radius. The results of the simulations are seen to converge to the analytic solution at later times, as more 
particles pass through the shock and the effects of numerical resolution become less important. 



At any radius inside the shock, and at all times, the ra- 
tio between the local dynamical time {3n/16Gpy^^ and the 
cooling time equals tv-s/isJTgKo. 

Equations 4.1, 4.2, 4.4, and 10 can now be solved to 
describe the post-shock flow once adequate boundary condi- 
tions at the center are imposed. As discussed by Bertschinger 
(1989), three different kinds of solutions can be identified 
according to the limiting behaviour of the solution near the 
center. The first type is a solution where the fiow stagnates 
at some finite radius. Infalling gas settles onto this surface, 
where the density formally diverges. In order to obey self- 
similarity the surface must move outwards at the same rate 
as the turnaround radius. This kind of solution thus requires 
a piston to move the surface outwards and is therefore of lit- 
tle physical applicability. 

The second type of inner solution corresponds to a flow 
that extends all the way to A = so that the local flow time, 
tfiow = r/v, near the center becomes much shorter than the 



cooling time. These solutions are referred to as "adiabatic" 
solutions, because cooling is unimportant for small A. The 
mass accretion rate near the center approaches a constant 
and the accretion speed diverges near the center. 

The third kind of solution, sometimes called the "eigen- 
solution", is the limiting adiabatic solution with minimum 
central mass accretion rate, or, equivalently, the limit of the 
family of stagnating solutions as the stagnation radius tends 
to zero. 

Each of these solutions is characterized by diflierent val- 
ues of the shock radius. As. Values of As similar to those 
obtained neglecting cooling correspond to solutions with 
stagnation points. As the shock radius moves inwards the 
stagnation point moves closer to the center and the solution 
transitions through the eigensolution to the adiabatic case. 

Figure 2 shows examples of these three different so- 
lutions for the case e = 1, Ko = 0.1. The solution with 
As = 0.23 (dotted line) has a stagnation point at Ao ~ 0.026 
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where the density diverges and the velocity becomes zero 
at the surface. There is no mass inside this radius, and the 
surface is pushed out by a piston to preserve similarity. As 
the shock radius is reduced to \s = 0.15 the flow extends 
all the way to the center and the infall velocity diverges 
there. The eigensolution corresponds to As « 0.1869. In this 
case the central velocity remains finite at the origin and the 
flow extends all the way to the center. In practice, we find 
this solution numerically by letting the stagnation point, Ao, 
approach zero. 

Noting that the solutions have approximately constant 
infall velocity in the inner regions, it is possible to determine 
the asymptotic slopes of the density and pressure profiles. 
As the velocity tends to a constant, D[\) —> DoX~^ and 
P(A) P()A^'^, where the normalising constants depend on 
the cooling parameter Kq. The asymptotic velocity is given 
by -7^0^^7(27 - 3). 

Figure 3 shows how the eigcnsolutions vary as a func- 
tion of the cooling efficiency parameter Kq. As the impor- 
tance of cooling increases the pressure support inside the 
shock decreases and the shock radius moves inwards. Per- 



haps counterintuitively, as cooling becomes more important 
the temperature inside the shock radius increases and the 
density decreases. This is because low entropy gas is "lost" 
to the central mass and, at fixed radius, low entropy gas 
is replaced by higher entropy gas that moves in from out- 
side. The "cooling flow" thus results in a net increase in gas 
entropy at a given radius. Figure 3 illustrates that cooling 
has a substantial effect on the structure of the system, and 
suggests that similarity solutions with cooling may provide 
a stringent test of the capabilities and accuracy of hydrody- 
namical codes. We pursue this issue next. 



3 COMPARISON WITH SPH SIMULATIONS 

As discussed in §1, the solutions derived in the previous 
section can be fruitfully confronted with the results of cos- 
mological hydrodynamical codes. This comparison is all the 
more interesting because the test case we discuss in the 
previous section captures many of the salient features of 
the galaxy formation process: gravitational collapse, pres- 
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e 




Kq 


As 


1 


8/9 


0.0 


0.3389 


1 


8/9 


0.1 


0.1858 


1 


8/9 


0.3 


0.0939 


2/3 


1 


0.0 


0.2899 


2/3 


1 


0.1 


0.1551 


2/3 


1 


0.3 


0.0733 


1/3 


4/3 


0.0 


0.1889 


1/3 


4/3 


0.1 


0.0996 


1/3 


4/3 


0.3 


0.0422 



Table 1. The dimensionless shock radius of the eigensolution 
for various choices of the shape parameter of the initial density 

perturbation, e, and of the dimensionless cooling coefficient, Xq- 
The time exponent of the turnaround radius, is also listed for 
each case. 



surization through shocks, radiative energy losses, cooling 
flows. Furthermore, because the solutions arc self similar in 
time, a single simulation can be examined at different times 
and convergence can be directly assessed. This is impor- 
tant becaxisc the importance of numerical resolution varies 
with time within a single sirrmlation. For example, the rmm- 
ber of particles within the shock radius increases with time, 
and the ratio between the smallest resolved radius and the 
turnaround radius decreases with time. Because the solution 
is unique, analyzing the deviations between analytic solution 
and numerical experiment at different times provides invalu- 
able insight into the role of numerical limitations and their 
consequence on the subsequent evolution of the system. 

We use the Smooth Particle Hydrodynamics code de- 
scribed by Navarro & White (1993), where details about the 
numerical procedure should be consulted. The initial setup is 
also similar to that described by these authors. We simulate 
a spherical region of an Einstein-de Sitter universe by lay- 
ing down 24, 257 particles homogeneously inside a sphere of 
radius R. Each particle is given an initial velocity consistent 
with unperturbed Hubble flow, and an external potential is 
added to mimic a point mass perturbation of mass equal to 
5% of the total mass of the sphere. The external potential 
is "softened" inside a fixed radius Rp = 0.1/? in order to 
prevent divergences, but is fully Kcplcrian outside Rp. This 
corresponds to the case e = 1 in eq.(l). 

The gravitational softening of each particle is also cho- 
sen to be equal to Rp. All particles have initially the same 
temperature, chosen to be much lower than the final virial 
temperature of the system in order to prevent hydrodynam- 
ical effects from becoming important before the gas turns 
around and passes through the shock. Two different sim- 
ulations were performed, one with Ko = and one with 
K{_) = 0.1. The sirrmlations are evolved until the turnaround 
radius encompasses half of the total number of particles. 
At the final time, about 25% of the particles have passed 
through the accretion shock. 

Figures 4 and 5 show the dimensionless profiles of den- 
sity, velocity, and temperature averaged in spherical bins 
of constant logarithmic width. The dimensionless mass en- 
closed inside each bin is also shown in the bottom right 



panel. Each panel shows the result of the simulations at 
three different times, corresponding to different numbers of 
particles within the shock radius: 1,838 (triangles), 3,648 
(squares), and 5,443 (circles) for the simulation without 
cooling and 1,655 (triangles), 3,371 (squares), and 4,734 
(circles) for the simulation with cooling. The analytic eigen- 
solutions corresponding to Ko =Q and Ko =0.1 are shown 
with a solid and dashed line, respectively. 

The simulation without cooling is in all respects similar 
to that reported by Navarro & White (1993), except for the 
fact that we use about ten times more particles. As shown in 
Figure 4, the density and velocity profiles agree remarkably 
well with the analytic solution inside the shock radius As. 
As expected, agreement with the analytic solution improves 
as more and more particles pass through the shock and the 
post-shock fiow becomes better resolved. This is especially 
noticeable in the mass panel, where it is seen that at early 
times (triangles) the simulation deviates from the analytic 
solution but that the system converges to the right solution 
at later times. Once about 5, 000 particles have gone through 
the shock the agreement between the solution and the ex- 
periment is remarkably good. The convergence towards the 
similarity solution is also convincingly demonstrated in the 
temperature panel, where the linear temperature scale ac- 
centuates the discrepancies near the center. 

One important conclusion from this analysis is that 
although numerical limitations compromise the results of 
the rmmcrical simulations at early times, these discrepan- 
cies have no major effect on the behaviour of the system 
at late times. The main deviations from the analytic solu- 
tions actually happen beyond the nominal shock radius. The 
shock is smoothed over several resolution lengths, and even 
at late times only outside approximately 2As the preshock 
fiow solution is recovered in the simulations. 

As discussed in §2.1 and illustrated directly by the 
dashed and solid lines in Figures 4 and 5, the similarity 
solution changes substantially when radiative cooling is in- 
cluded. The shock radius moves inwards, the post-shock 
temperature increases, the infall velocity is non zero all the 
way to the center, and the density decreases at all radii be- 
cause of the accumulation of cold material at the center. 

Figure 5 shows that the numerical results match very 
well these predicted changes in the post-shock region. In- 
side As = 0.1872 the temperature, mass, density, and veloc- 
ity profiles are almost indistinguishable from the similarity 
solutions. Remarkably, this is the case even when only 7% 
(« 1,600 particles) of the mass of the system has passed 
through the shock, and the agreement is seen to improve 
as more and more particles go through the shock. As noted 
above, the main shortcoming of the simulation regards the 
width of the shocked region: particles are seen to respond 
to the shock as far out as ~ 3 times the nominal shock 
radius. It is encouraging, however, to note that the overall 
trend is correct, and that the volume of the shocked region 
is substantially smaller than in the case where cooling is 
neglected. 



4 DISCUSSION 

One important goal of the comparison between simulations 
and similarity solutions presented above is to assess the reli- 
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ability of numerical techniques currently being used to sim- 
ulate the formation of galaxy-sized structures in the uni- 
verse. The present study indicates that, to a large degree, 
SPH methods give results that are consistent with analytic 
solutions in test cases that involve some of the major in- 
gredients believed to play a significant role during galaxy 
formation. In particular, our results show that SPH codes 
can reproduce faithfully the changes in temperature, density, 
and velocity that are associated with strong "cooling flows" 
onto a central perturbation. The mass inflow rates are also 
accurately reproduced, and there is no indication that nu- 
merical limitations lead to an undue increase or decrement 
in the amount of mass that cools and flows to the central 
object. 

This is important because it has been argued that SPH- 
like treatments of numerical hydrodynamics may cause an 
artiflcial "overcooling instability" that exaggerates the im- 
portance of radiative cooling losses and leads to the forma- 
tion of excessively massive gas concentrations at the center 
of dark halos (Thacker et al 1998) . Our study shows that this 
instability is not present in the test cases we present here. 
The rate at which gas cools and gets accreted onto a central 
clump is proportional to the mass accretion rate through the 
shock radius and is consistent with that expected from the 
similarity solution. 

There are, however, important differences between the 
cases considered in our study and in that of Thacker et al. 
These authors consider gas cooling from a hot gaseous halo 
in hydrostatic equilibrium within a dark halo and report a 
strong dependence of the cooled mass on the numerical res- 
olution of SPH simulations. They also consider a more re- 
alistic cooling function that varies with temperature and 
density in very different ways than the self-similar cooling 
function we adopt here. Because of these caveats, it may be 
premature to argue either for or against their flndings on the 
basis of the simulations presented here. It should be possible, 
however, to test their arguments using the cooling wave sim- 
ilarity solutions derived by Bertschinger (1989). These are a 
much closer analogue to the case considered by Thacker et al 
and direct comparison between numerical experiments and 
Bertschinger's analytic solutions should provide a definitive 
assessment regarding the effects of numerical resolution on 
the cooling and accretion of gas at the center of dark halos. 
We plan to carry out this comparison in the near future. 



5 SUMMARY 

We have derived similarity solutions that describe the spher- 
ical collapse, shock, and settling of collisional gas from scale- 
free pPI-tTirhatinng in an Kinatpin-Hp Sittgr irniwrap Onr 



tree pprturhatmng in an H.inatPin-rlP Nittpr nniwrap j 

study ekterids prior work on the subject by taking into 
count t he full cffccta of energy lo33 due to radiative cooling 



process es, under the simplif^dng assumption that the cooling 
function is a simple power-law of density and temperature. 
This choice ensures that the time evolution of the system is 
self-similar by requiring that the cooling time of the system 
at all times is a fixed fraction of the age of the universe. 
The solutions take into account many of the processes that 
are relevant to the assembly of the baryonic component of 
galaxies in a cosmological scenario: gravitational infall, en- 
ergy dissipation through shocks, pressure gradients, radia- 



tive energy losses, and cooling flows. Analytic solutions such 
as the ones outlined here are invaluable to assess the reliabil- 
ity and diagnose the shortcomings of numerical techniques 
currently being used in cosmological simulations. 

The tests we present here show that SPH simulations 
reproduce the analytic solutions very well. No substantial 
deviations from the predicted central mass accretion rates 
or from the temperature, density, and velocity profiles are 
observed inside the shock radius when cooling is included. 
The region affected by the shock is, however, is much larger 
than predicted: effects from the shock are seen as far out as 
2 or 3 times the nominal shock radius. Although this does 
not seem to affect adversely the post-shock behaviour of the 
gas under the simplifying conditions we adopt, it may have 
unwanted consequences in cases with more complex infall 
geometry (such as mergers) or that involve a cooling func- 
tion with a more sensitive dependence on temperature and 
density than assumed here. We hope that the work reported 
here will encourage further efforts to test and improve the 
numerical treatment of the hydrodynamics of galaxy forma- 
tion. 
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